/******************************************************************************
 * mySCPsolver.cpp: Implementation for My Steiner Constrained Path Solver class.
 *
 * Author: Carlos Eduardo de Andrade <andrade@ic.unicamp.br>
 *
 * (c) Copyright 2011 Institute of Computing, University of Campinas.
 *     All Rights Reserved.
 *
 *  Created on : Oct 01, 2011 by andrade
 *  Last update: Oct 02, 2011 by andrade
 *
 * This software is licensed under the Common Public License. Please see
 * accompanying file for terms.
 *****************************************************************************/

#include <iostream>
using namespace std;

#include <lemon/list_graph.h>
using namespace lemon;

#include <typeinfo>
#include "mySCPsolver.hpp"
#include "../timer.hpp"
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sstream>
#include "gurobi_c++.h"

//-------------------------[ Default Construtor ]----------------------------//

mySCPsolver::mySCPsolver():
		SteinerConstrainedPathSolver(),
		ze() {}

//-------------------------[ Default Destructor ]----------------------------//

mySCPsolver::~mySCPsolver() {}

//-----------------------------[ Set Instance ]------------------------------//

// Se for sobrecarregar esta funcao, tome cuidado de limpar o grafo antigo.
// Recomendo chamar o mehtodo da classe base primeiro.
bool mySCPsolver::loadInstance(const char* filename, double delayParam) {
    return SteinerConstrainedPathSolver::loadInstance(filename, delayParam);
}

void mySCPsolver::cortafolhas(int folha, ListGraph &graph){
	while(folha == 0){
		folha = 1;
		for (int id = 0; id <= SteinerConstrainedPathSolver::graph.maxNodeId(); ++id) {
			for (ListGraph::NodeIt n(SteinerConstrainedPathSolver::graph); n != INVALID; ++n) {
			    int cnt = 0;

			    //Conta o numero de arestas incidentes no nó n
			    for (ListGraph::IncEdgeIt e(SteinerConstrainedPathSolver::graph, n); e != INVALID; ++e) {
			      cnt++;
			    }
				if (cnt == 1 && (SteinerConstrainedPathSolver::terminal[n] == false) && (SteinerConstrainedPathSolver::extremos[n] == false) ){
			    	SteinerConstrainedPathSolver::graph.erase( n );
			    	folha = 0;
			    	std::cout << "deg(" << SteinerConstrainedPathSolver::graph.id(n) << ") = " << cnt << std::endl;
			    }
			  }
		}
	}
}

string itos(int i) {stringstream s; s << i; return s.str(); }

void allsubsets(int s[],int p,int k,int t[],GRBVar **vars,GRBModel &model,int q=0,int r=0) {
    if(q==k)
    {
    	GRBLinExpr expr = 0;
        for(int i=0;i<k;i++) {
        	        	
        	for (int j=0;j<k;j++) {

        		expr += vars[ t[i] ][ t[j] ];
        	}
        }
        model.addConstr(expr <= k - 1, "subtour_"+itos(k));
    }
    else
    {
        for(int i=r;i<p;i++)
        {
            t[q]=s[i];
            allsubsets(s,p,k,t,vars,model,q+1,i+1);
        }
    }
}

int proximoNo(int noID,int n,int final, GRBVar **vars){

	if (noID == final){
		return -1;
	}
	else {
	
		for (int i = 0; i < n; i++) {
			if (vars[noID][i].get(GRB_DoubleAttr_X)) {
				return i;
			}
		}

		for (int i = 0; i < n; i++) {
			if (vars[i][noID].get(GRB_DoubleAttr_X)) {
				return i;
			}
		}	
	}	

}

//--------------------------------[ Solve ]----------------------------------//

mySCPsolver::ResultType mySCPsolver::solve(const double max_time) {
	//Resolva aqui de forma exata

	Timer timer;
    timer.restart();
	while(timer.elapsed() < max_time - 1){

		ListGraph::Node aux;
		ListGraph::Node aux2;
		ListGraph::Node s; // nó inicial
		ListGraph::Node t; // nó final

		int sIsDefined = 0;
		int sID,tID,result;
		int found_any = 0;
		std::vector<ListGraph::Node> current_solution;
		double value_current_sol = 0.0;

	    GRBEnv env = GRBEnv();
	    GRBModel model = GRBModel(env);

	    int n = countNodes (graph);
	    GRBVar **vars = new GRBVar*[n]; //matriz de adjacencia das arestas
	    GRBLinExpr steinerConstr1[n];
	    GRBLinExpr steinerConstr2[n];
	    double delayMatrix[n][n]; // matriz auxiliar que conterá os delays entre cada nó

		std::cout << "Número de nós:" << std::endl;
	 	std::cout << countNodes (graph) << std::endl;

	 	//Adiciona todos as arestas do grafo como variáveis binárias
	 	//1 se o nó faz parte da solução
	 	//0 caso contrário
	 	for (int id = 0; id <= SteinerConstrainedPathSolver::graph.maxNodeId(); ++id) {
	 		vars[id] = model.addVars(countNodes (graph));
	 	}
	 	model.update();

	 	for (int i = 0; i < n; i++) { // inicialização da matriz de adjacência
	 		for (int j = 0; j < n; j++) {
	 			aux = SteinerConstrainedPathSolver::graph.nodeFromId(i);
	 			aux2 = SteinerConstrainedPathSolver::graph.nodeFromId(j);

	 			ListGraph::Edge e = findEdge(graph,aux,aux2);
	 			if(e != lemon::INVALID) {
	 				vars[i][j].set(GRB_DoubleAttr_Obj, length[e]);
	 				delayMatrix[i][j] = delay[e];
	 			}
	 			else {
	 				vars[i][j].set(GRB_DoubleAttr_Obj, 99999); //pog pra setar custo infinito caso a aresta nao exista
	 				delayMatrix[i][j] = 99999;
	 			}

	 			vars[i][j].set(GRB_CharAttr_VType, GRB_BINARY);
	 			vars[i][j].set(GRB_StringAttr_VarName, "x_"+itos(i)+"_"+itos(j));
	 		}
	 	}
	 	model.update();

	 	for (int i = 0; i < n; i++) { //nesse loop fixa o i e varia o j
			aux = SteinerConstrainedPathSolver::graph.nodeFromId(i);
	 		
	 		if (SteinerConstrainedPathSolver::extremos[aux]==true) {
	 			
	 			if (!sIsDefined) {//SUM(Xsj) = 1
	 				s = aux;
	 				sID = i;
	 				sIsDefined = 1;

	 				GRBLinExpr expr = 0;

	 				for (int j = 0; j < n; j++)
	 					expr += vars[i][j]; // somatória de todos as arestas saindo do vértice inicial
	 				model.addConstr(expr == 1, "deg1_ex_sj_"+itos(i));
	 			} 
	 			else { //SUM(Xtj) = 0
	 				t = aux;
	 				tID = i;

	 				GRBLinExpr expr = 0;

	 				for (int j = 0; j < n; j++)
	 					expr += vars[i][j]; //somatória de todos as arestas saindo do vértice final
	 				model.addConstr(expr == 0, "deg1_ex_tj_"+itos(i));
	 			}
		 	}

		 	else if (SteinerConstrainedPathSolver::terminal[aux]==true){
				// Nós terminais k devem ter grau 1
				// SUM(Xkj) = 1
		 		GRBLinExpr expr = 0;

		 		for (int j = 0; j < n; j++)
		 			expr += vars[i][j];
		 		model.addConstr(expr == 1, "deg2_ex_kj_"+itos(i));	 		
		 	}

		 	else{ //Vértice de steiner l
		 		//SUM (Xli) <= 1
				GRBLinExpr steiner2 = 0;
		 		for (int j = 0; j < n; j++)
		 			steiner2 += vars[i][j]; //Somatoria de todas as arestas saindo de um vertice de steiner
		 		model.addConstr(steiner2 <= 1, "deg2_ex_lj_"+itos(i));
		 		steinerConstr2[i] = steiner2;
		 	}
		}
		model.update();

		
		for (int j = 0; j < n; j++) { // nesse loop fixa o j e varia o i

			aux = SteinerConstrainedPathSolver::graph.nodeFromId(j);

			if (SteinerConstrainedPathSolver::extremos[aux]==true) {

				if (aux == s) {
		 			// SUM(Xjs) = 0
		 			GRBLinExpr expr = 0;

			 		for (int i = 0; i < n; i++)
			 			expr += vars[i][j];
			 		model.addConstr(expr == 0, "deg1_ex_js_"+itos(j));
		 		}
		 		else if (aux == t) {
		 			// SUM(Xjt) = 1

		 			GRBLinExpr expr = 0;
		 			for (int i = 0; i < n; i++)
		 				expr += vars[i][j];
		 			model.addConstr(expr == 1, "deg1_ex_jt_"+itos(j));
		 		}
		 	}
			else if (SteinerConstrainedPathSolver::terminal[aux]==true){
				// SUM (Xjk) = 1
		 		GRBLinExpr expr = 0;

		 		for (int i = 0; i < n; i++)
		 			expr += vars[i][j];
		 		model.addConstr(expr == 1, "deg2_ex_jk_"+itos(j));	
			}

			else{ //Vértice de steiner l
		 		//SUM (Xjl) <= 1
				GRBLinExpr steiner1 = 0;
		 		for (int i = 0; i < n; i++)
		 			steiner1 += vars[i][j]; //Somatoria de todas as arestas chegando num vertice de steiner
		 		model.addConstr(steiner1 <= 1, "deg2_ex_jl_"+itos(j));
		 		steinerConstr1[j] = steiner1;

		 		GRBLinExpr expr = steinerConstr1[j] - steinerConstr2[j];
		 		model.addConstr(expr == 0, "diff0_ex_jl_"+itos(j));
		 	}
		}
		model.update();

		//Restrição do delay máximo
		GRBLinExpr expr = 0;
		for (int i = 0; i < n; i++) { 
			for (int j = 0; j < n; j++) {
				expr += vars[i][j] * delayMatrix[i][j];
			}
		}
		model.addConstr(expr <= SteinerConstrainedPathSolver::delayLimit, "max_delay");
		model.update();

		//Restrições de subotour
		int v[n];
		for (int i = 0; i < n; i++) {
			v[i] = i;
		}

		for (int i = 2; i < n; i++) {
			int v2[i];
			allsubsets(v,n,i,v2,vars,model); 
		}

		//Tenta achar alguma solução factível pelas heurísticas
		mySCPsolver::solveFast(0.2 * max_time);

		if (best_solution.size()>0) { 
			//ACHOU solução pela heurística
			//Seta novos valores iniciais para o gurobi

			double hsolution = 0.0; //custo da solução obtida pela heurística

			for(unsigned i = 0; i < best_solution.size()-1; ++i){
				int nsrc = graph.id(best_solution[i]);
				int ndes = graph.id(best_solution[i+1]);

				// Seta os valores das arestas como 1 caso elas também façam parte da solução da heurística,
				// Caso contrário seta para 0
				vars[ nsrc ][ ndes ].set(GRB_DoubleAttr_Start,1); 
				vars[ ndes ][ nsrc ].set(GRB_DoubleAttr_Start,1);

				ListGraph::Edge e = findEdge(graph,best_solution[i],best_solution[i+1]);
				if(e != lemon::INVALID)
					hsolution += length[e]; //soma o custo da solução obtida pela heuristica
				}

				model.getEnv().set(GRB_DoubleParam_Cutoff, hsolution); //seta como ub inicial para o gurobi
			

			for (int i = 0; i < n; i++) { 
				for (int j = 0; j < n; j++) {
					
					if( vars[ i ][ j ].get(GRB_DoubleAttr_Start) !=1  )
						vars[ i ][ j ].set(GRB_DoubleAttr_Start,0);

				}
			}
		}


		model.update(); model.getEnv().set(GRB_IntParam_OutputFlag, 0);

		model.optimize();

		//Se achou solução pelo gurobi
		if (model.get(GRB_IntAttr_SolCount) > 0) {
			found_any = 1;
			//std::cout << "UHUUUUUL!" << std::endl;
		
			//Constroi o vetor de solução
			int inicio = sID;
			int prox;
			current_solution.push_back(s);
			while (prox != -1) {
				prox = proximoNo(inicio,n,tID,vars);
				if (prox != -1) {
					aux = SteinerConstrainedPathSolver::graph.nodeFromId(prox);
					current_solution.push_back(aux);
					inicio = prox;
				}
			}

			value_current_sol = 0;
			for(unsigned i = 0; i < current_solution.size()-1; ++i) {
				ListGraph::Edge e = findEdge(graph,current_solution[i],current_solution[i+1]);
				if(e != lemon::INVALID)
					value_current_sol += length[e];
			}

			result = checkSolution(current_solution, value_current_sol);
			printf("EXACT RESULT:%d\n",result);

			if(result == OK && value_current_sol < solution_value) {
				solution_value = value_current_sol;
				best_solution.clear();
				best_solution.resize(current_solution.size());
				for(unsigned i = 0; i < current_solution.size(); ++i){
					best_solution[i] = current_solution[i];
				}
				found_any = 1;
				printf("Solução atual:%f\n",value_current_sol);
			}
		}


		if(found_any)return mySCPsolver::EXACT_SOLUTION;
	    return mySCPsolver::EXACT_NO_SOLUTION;

	}
}

//--------------------------------[ Solve ]----------------------------------//

mySCPsolver::ResultType mySCPsolver::solveFast(const double max_time) {

    Timer timer;
    timer.restart();

    //Solucao atual
    std::vector<ListGraph::Node> current_solution;
    double value_current_sol = 0.0;
    int numNodeSol, result = 1, found_any = 0;
	int vez = 0;
	int x = 0;

	//random seed
	//srand(time(NULL));

	//Enquanto houver tempo
	while(timer.elapsed() < max_time - 1){
		result = 1;

		// Verifica se existe algum nó de grau 1 que não seja terminal
		// Se houver ele é removido do grafo original
		int folha = 0;
		if (vez == 0)
			cortafolhas(folha,graph);
		vez = 1;

		if (x == 0) {
			algoritmoConstr alg1;
			value_current_sol = alg1.construirSolucao(current_solution);
			result = checkSolution(current_solution, value_current_sol);
			printf("RESULT:%d\n",result);
			x = 1;
		}
	
		if(result == OK && value_current_sol < solution_value) {
			solution_value = value_current_sol;
			best_solution.clear();
			best_solution.resize(current_solution.size());
			for(unsigned i = 0; i < current_solution.size(); ++i){
				best_solution[i] = current_solution[i];
			}
			found_any = 1;
			printf("Solução atual:%f\n",value_current_sol);
		}

		//HEURISTICA de trocar a ordem dos vértices para diminuir a distancia final
		//Só pode rodar se a entrada for alguma solução factível
		if (result == OK || result == INCORRECT_DELAY) {
			algoritmoBuscaLocal alg;
			double solucao = alg.buscaLocal(current_solution, graph, value_current_sol);		
			
			if (solucao == 1) { //só entra aqui se a solução foi melhorada
				value_current_sol = 0;
				for(unsigned i = 0; i < current_solution.size()-1; ++i) {
					ListGraph::Edge e = findEdge(graph,current_solution[i],current_solution[i+1]);
					if(e != lemon::INVALID)
						value_current_sol += length[e];
				}

				result = checkSolution(current_solution, value_current_sol);
				if (solucao != 0 && result != OK) {
					//printf("OUtro valor:%f\n",value_current_sol);
					//printf("Result:%d\n",result);
				}

				if(result == OK && value_current_sol < solution_value) {
				//if(value_current_sol < solution_value) {
					solution_value = value_current_sol;
					best_solution.clear();
					best_solution.resize(current_solution.size());
					for(unsigned i = 0; i < current_solution.size(); ++i){
						best_solution[i] = current_solution[i];
					}
					found_any = 1;
					printf("Solução atual melhorada:%f\n",value_current_sol);
				}
			}
		}

		//cout << "\n";
		//vamos procurar a proxima solucao
		current_solution.clear();
		value_current_sol = 0;
	}
	ub = solution_value;
	if(found_any)return mySCPsolver::FAST_HEURISTIC_SOLUTION;
    return mySCPsolver::FAST_HEURISTIC_NO_SOLUTION;
}
